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Quantum Monte Carlo simulation of three-dimensional Bose-Fermi mixtures 
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Exploratory simulations of Bose-Fermi mixtures on the three-dimensional optical lattice at finite 
temperature are performed by adopting the lattice quantum chromodynamics technique. We analyze 
the bosonic superfluid transition and its dependence on the strength of the boson-fermion coupling. 
The particle densities and the pair occupancies are also studied to understand the effect of the boson- 
fermion coupling to the microscopic properties of the system. Effect of the induced fermion-fermion 
interaction by the boson density fluctuation is clearly seen. 
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I. INTRODUCTION 

Many-body mixture of bosons and fcrmions is an in- 
teresting quantum system. It is not only useful to study 
the interplay between different quantum statistics in con- 
densed matter physics and in atomic physics but 
also important for understanding the hadron-quark phase 
transition in nuclear and particle physics 0. Recently, 
many kinds of Bose-Fermi mixtures have been realized 
in ultracold atomic experiments. In particular, the Bose- 
Fermi mixtures such as ^^Rb-^OK i^OYb-^^^yb and 
^^"'Yb-^^'^Yb @ can be trapped in three-dimensional op- 
tical lattices. 

Despite of its importance, ab initio quantum Monte 
Carlo calculations of the Bose-Fermi Hubbard model, 
which is a model system of the Bose-Fermi mixture on the 
optical lattice, have been limited to one spatial dimension 
[7l-[l3j. (See also a recent review These simulations 

were done with improved versions of the world-line for- 
malism [Tsl - flTj . The world-line formalism is an exact 
scheme in one spatial dimension while it has the fermion 
sign problem in higher dimensions. 

In this paper, we apply the method of lattice Monte 
Carlo simulation of quantum chromodynamics (QCD) 
[isj to overcome the fermion sign problem in the Bose- 
Fermi mixture of ultracold atoms on an optical lattice in 
three spatial dimensions at finite temperature. The lat- 
tice QCD simulation is an established numerical frame- 
work for the strongly coupled regime of QCD; the mix- 
ture of gluons (bosons) and quarks (fermions) can be 
treated in any dimension without fermion sign problem as 
long as the fermion determinant is real and semi-positive. 
(For relativistic fermions with finite chemical potential, 
the fermion determinant becomes complex, while the 



chemical potential for non-relativistic fermions does not 
cause the problem.) 

This paper is organized as follows. In Sec. II, we in- 
troduce the model action and discuss the basic setup of 
our simulations. In Sec. Ill, we analyze the bosonic su- 
perfluid transition and its dependence on the strength of 
the boson-fermion coupling. The particle densities and 
the pair occupancies are also studied to understand the 
effect of the boson-fermion coupling. In Sec. IV, we sum- 
marize this study. A part of the preliminary results was 
reported in Ref. [l9| . 



II. FORMULATION 

We consider one-component boson field $(2;, t) and 
two-component fermion field 5'-|-(af, r) and '^^{x^t). The 
partition function Z is then written as a path integral in 
terms of the Euclidean action S: 



det A> dot A'4,e-^« . (1) 



In the second line, the bilinear form of the fermion fields 
are integrated out, so that Z becomes a functional only 
of the boson field. 

We adopt the single-band Bose-Fermi Hubbard model 
whose lattice action reads [l9l-[2ll 



S — Sb ^ Sp -\- S 



BF, 



(2) 



Sb = ['^*{x,t){<^>{x,t) - e^«$(f,T - 1)} -^ii3{$*(x,T)$(x + e;-,T) + $*(x,t)$(x - e;-,T)} 

X.T j = l 

+C/Bs$*(a?,r)$(f,T){$*(a-,T)$(x,T) - 1}] (3) 
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Sf = J2 [KiS,r){'i',ix,T)-e''-'i'^{x,T-l)}-Y,tF{KiS,r)'i',ix + ej,T) + ^;ix,T)^ 



SbF = ^BF<f*{x,T)<i>{x,T)^l{x,T)^^{x,T ^l) 



(5) 



r 



Here, /is (mf) is the boson (fermion) chemical poten- 
tial, (tp) is the spatial hopping parameters for the 
boson (fermion). Also, Ubb and Ubf are the strength of 
boson-boson interaction and boson-fermion interaction, 
respectively. We do not consider the fermion-fermion in- 
teraction for simplicity, i.e. Uff = 0. The temporal 
hopping structures of the chemical potential and the in- 
teraction are based on Rcfs. (20l. [2l|. 

Note that we have absorbed the spatial lattice spacing 
Gs and the temporal lattice spacing ar into the definition 
of the fields and the coupling parameters. Here as should 
be identified with a finite lattice constant of an optical 
lattice, while Or should be eventually taken to be zero 
to recover the continuum limit in the temporal direction. 
In the following discussion, all dimensional quantities are 
scaled by the lattice constant. 

Our lattice action has the boson sign problem which 
appears in non-relativistic systems. The boson action 
is complex in general because the temporal hopping 
term is given as t){$(x, r) — $(.t, r — 1)} = 

Nr E,.{l-e"'"''^''^}**(^, iy)^{x, v) with the boson Mat- 
subara frequency ujb{i^) = ^vnT. To avoid this boson 
sign problem, we adopted the zero-frequency approxima- 
tion, in which the boson field is projected onto the zero- 
frequency mode as r) = $(a?, = 0) = This 
approximation is validated if the system is either in the 
critical region or in the high-temperature limit, where the 
static mode of the boson field, v = gives a dominant 
contribution. Furthermore, effect of the non-zero Mat- 
subara modes can be estimated by using the reweigliting 
method [2^ . 

We can prove that the fermion determinant is real 
and positive. Let us consider the fermion action in 
the frequency space. The matrix K„ becomes diago- 
nal in the frequency space. The complex factor is only 
^-iu}f{v) fermion Matsubara frequency uif{v) = 

{2v — 1)t:T. The fermion determinant is 



B) 



= Y\_ [{A + Bcosujf{v)Y + B"^ sm^ ijJf{v)] 

> , (6) 

where A and B are real numbers. Therefore, there arises 
no sign problem originating from the fermion determi- 
nant. (This is in contrast to the Bose-Fermi system with 
a Yukawa coupling discussed e.g. in Ref. [23|.) Even if 
the fermion-fermion interaction is introduced, there is no 
fermion sign problem: After the Hubbard- Stratonovich 



transformation of the fermion-fermion interaction, K^j is 
still a real matrix with auxiliary fields, so that the two- 
component fermion determinant becomes semi-positive, 
det i^t det = (det K-^f > 0. 

Our numerical simulation procedure is based on the 
hybrid Monte Carlo algorithm with the pseudofermions 
2J] , which is widely used in modern simulations of lattice 
field theories. For preconditioning the fermion matrix, 
we perform a Fourier transformation from the imaginary 
time T to the fermion Matsubara frequency llJf{i')- In 
this study, we fixed the hopping parameters Ib = tp = 
0.01 and the repulsive boson-boson interaction Ubb = 
0.1 for the first trial. We performed the simulation in 
a spatial lattice volume V = = 10^ with periodic 
boundary conditions. We varied temperature T = 1/A'^^ 
by changing the temporal lattice size as Nr = 20, 30, 40, 
60, and 100. 

In the grand-canonical approach, particle number den- 
sities are related to chemical potentials and other param- 
eters. Therefore, we monitor the boson number density 
ub and fermion number density np in our simulations: 

UB = ($*(x)$(x)), (7) 
np = {mx,T)^>^{x,r)) = (^l{x,T)^^ix,r)). (8) 

Note that, in the mean-field level, the boson-fermion cou- 
pling Ubf induces effective chemical potentials: 



Ml? 



^J■B - UBpnp, 
fip - UBpnB- 



(9) 
(10) 



This implies that attractive boson-fermion interaction 
{Ubf < 0) increases the effective chemical potentials and 
hence the particle number densities. Our numerical sim- 
ulations given below indeed show such behavior. 



III. NUMERICAL RESULTS 

To see the effect oi Ubf on the Bosc-Einstein conden- 
sation, let us first examine the boson propagator as a 
function of the distance r ~ \x ^ y\, 



(11) 



The boson propagator Gb(?') with Ubf = and fip ^ 
jjLB = is shown in Fig. [1] Although the statistical er- 
rors are shown by the error-bars in the figure, they are 
smaller than the size of the symbols and practically in- 
visible. At high temperature (T > 0.02), G'B(r) drops to 
zero exponentially, while at low temperature (T < 0.02), 
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Gsir) approaches to a constant value for large r, i.e. a 
signature of the the off-diagonal long-range order. 

In Fig. [2j we draw the boson condensate fraction de- 
fined by 



Rl 



GB{r = Ns/2) 



(12) 



for different values of T and Ubf at np = hb = 0. The 
blue region corresponds to the superfluid phase. The 
figure shows that the superfluid region is enlarged as 
the boson-fermion attraction {Ubf < 0) becomes strong. 
This is because /i^ and hence the boson number density 
ub increase associated with the increase of the boson- 
fermion attraction. 

To see the effect oi Ubf on particle number profiles in 
more details, let us introduce normalized pair occupan- 
cies which is a measure to find a pair at the same lattice 





FIG. 3: The particle number densities and pair occupancies 
at fiB = hf ~ and T — 0.05 as a function of Ubf- 



FIG. 1: The boson propagator GB{r) as a function of r = 
|a; — j7| for different values of temperature T with Ubf ~ 
and fiF = fJ.B ~ 0. 
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FIG. 2: The bosonic condensate fraction Rb in the T-Ubf 
plane at = Ms = 0. 



site; 



Dbf 
Df 
Db 



{nB{x)nFa{x)) 
{nB{x)){nFa{x)) 

{nF^{x)nFi{x)) 
{nFtix)}{nFiix)) 

1 {nB{x)nB{x)) 

2 {nB{x))^ " 



(13) 
(14) 
(15) 



Here the number density operators are defined as 
nB{x) = $*(f)$(x) and npaix) = 4"* (f , r)^'cr(f , r). 
The coefficient I /2 in Z?b is a symmetry factor for iden- 
tical particles. The pair occupancies are experimentally 
observed on a three-dimensional optical lattice @ . 

In Fig. [3l we show the particle number densities (ns 
and uf) and the pair occupancies {Dbf, Dp, and Z3_b) as 
a function of Ubf- The chemical potentials are fixed at 
Ms = M-F = as in the previous figures, while the temper- 
ature is chosen to be T = 0.05 (upper boundary in Fig.[2|). 
As the boson-fermion attraction becomes stronger, both 
Ub and np increase, which is consistent with the behav- 
ior expected from the effective chemical potentials as we 
have already mentioned. Also, the pair occupancy of dif- 
ferent (same) species tends to be enhanced (suppressed) 
as Ubf increases, which is naturally expected from the 
contact nature of the boson-fermion interaction. Note 
that Db is always below unity, since we have taken re- 
pulsive boson-boson interaction. 

In Figs, m and [5j we show the particle profiles as a 
function of the fermion chemical potential fj,F for char- 
acteristic boson-fermion interactions; an attractive case, 
Ubf = —0.05, and a repulsive case, Ubf = +0.05. Other 
parameters are the same as those adopted in Fig. [31 By 
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FIG. 4: The particle profiles as a function of /if for attractive 
boson- fermion interaction Ubf = —0.05. Other parameters 
are the same with Fig. [S] 



pair occupancy Dbf is reduced (enhanced) from unity 
in the repulsive (attractive) case as naturally expected 
from the contact nature of the boson-fermion interaction. 
The fermion number density approaches to unity, when 
fi'^p in Eq. ([TO)) becomes large either by large attractive 
boson-fermion interaction or by large fermion chemical 
potential. If fermions fully occupy all the lattice sites, 
both Dbf and Df become unity, since the numerators 
of Eqs. dTS]) and factorize. 

Although the bare fermion-fermion interaction is 
switched off in our present simulations (Uff = 0), there 
arises induced fermion-fermion attraction due to the cou- 
pling of the fermions to the boson density fluctuation 
with the mag nitudeof 0(C/|p,/C/bb) Q- Because of this, 
the fermion pair occupancy Df is enhanced irrespective 
of the sign oi Ubf as seen in Figs. S] and [5] Also, the 
boson pair occupancy Db changes slightly because the 
boson number density ub changes. The boson pair oc- 
cupancy Db approaches to unity as the boson becomes 
dilute. 




FIG. 5: The particle profiles as a function of /if for repulsive 
boson-fermion interaction Ubf = -f0.05. Other parameters 
are the same with Fig. [S] 



varying fiF in the interval [—0.2, -1-0.2], the fermion num- 
ber density is changed from zero (i.e. complete va- 
cancy) to unity (i.e. full occupancy). The boson num- 
ber density ns is an increasing (decreasing) function of 
Hf in the attractive (repulsive) case, which is consis- 
tent with the behavior expected from the effective chem- 
ical potential in Eq. Also, the boson-fermion 



IV. SUMMARY 



In summary, wc have performed a quantum Monte 
Carlo simulation of Bose-Fermi mixture at finite tem- 
perature on the three-dimensional optical lattice. The 
sign problem from the boson hopping term is avoided by 
taking only the lowest Matsubara mode, which is valid 
near the critical point or at high temperature. We have 
analyzed the off-diagonal long-range order of the boson, 
the particle number densities, and pair occupancies for 
attractive and repulsive boson-fermion interaction. We 
found that the region of the Bosc-Einstein condensa- 
tion is extended as the boson-fermion attraction becomes 
strong. The particle number profiles obtained from the 
simulations can be understood in terms of the effective 
chemical potentials and induced fermion-fermion interac- 
tion driven by the boson-fermion interaction. 

The present exploratory study is a starting point of 
more expanded simulations with the ingredients includ- 
ing (i) the bare fermion-fermion interaction, (ii) finite size 
scaling analysis with several different lattice volumes, (iii) 
continuum limit in the temporal direction, (iv) taking 
into account the non-zero Matsubara mode for bosons, 
and (v) mapping of the phase diagram to the fixed num- 
ber of bosons and fermions. From the physics point of 
view, there are various interesting future problems to 
be examined with our approach, e.g. the superfluidity 
of fermion pairs induced by the density fluctuation of 
the background bosons with the boson-fermion coupling, 
and the properties of the boson-fermion mixture in the 
strongly coupled regime. 
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